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Abstract: Monte Carlo simulation of a 2+1 dimensional model of voltage-biased bi¬ 
layer graphene, consisting of relativistic fermions with chemical potential fi coupled to 
charged excitations with opposite sign on each layer, has exposed non-canonical scaling 
of bulk observables near a quantum critical point found at strong coupling. We present 
a calculation of the quasiparticle dispersion relation E{k) as a function of exciton source 
j in the same system, employing partially twisted boundary conditions to boost the 
number of available momentum modes. The Fermi momentum kp and superfluid gap A 
are extracted in the j ^ 0 limit for three different values of p, and support a strongly- 
interacting scenario at the Fermi surface with A ~ 0{p). We propose an explanation 
for the observation p, < kp m terms of a dynamical critical exponent z <1. 
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1 Introduction 


There are very few many-body systems permitting Monte Carlo simulation without the 
need to confront a Sign Problem. In Ref. [1] we introduced a new member to this class, 
based on an effective theory of bilayer graphene in which charge carrying excitations 
are modelled as Nf = 4 relativistic fermions moving in a 2d plane. The introduction of 
chemical potential /x is via a bias voltage in the perpendicular direction which induces 
equal densities of electrons on one layer and holes on the other; this is analogous to isospin 
chemical potential in QCD and yields a real positive fermion determinant amenable to 
orthodox Monte Carlo methods. The simulation was performed in the vicinity of a 
quantum critical point (QCP) found at strong coupling, and the main result was the 
demonstration that the ground state is a superfluid formed by condensation of electron- 
hole exciton pairs, and that the response to chemical potential is governed by the non- 
canonical scaling forms ([7]) and (E]) given in Sec. El 

The model [1] as originally devised for bilayer graphene |2] is artificial in a few 
respects. Firstly, the description in terms of iVj = 4 relativistic species (ie. Nf = 2 elec¬ 
trons and Nf = 2 holes) is only justihed by the band structure of the tight-binding model 
in the presence of an inter-layer “skew” coupling breaking the trigonal symmetry of the 
underlying lattice [3] • Secondly, the interaction between charge densities is simplihed to 
be a local four-fermi contact, although a more realistic unscreened Coulomb interaction 
can be modelled with the introduction of a third spatial lattice direction to capture the 
electrodynamics |1]. Finally, intra- and inter-layer interactions have the same coupling 
strength, as a necessary condition of keeping the fermion determinant real. Nonetheless, 
it shares the essential features of a more general model for double-layer graphene systems 
in which there is some hybridization permitting interlayer tunnelling [5]. With Nf = 1 
this approach is also applicable to surface states of topological insulators, motivating 
study with variable Nj [6]. 

The results of [1] were interpreted in terms of strong interaction effects at a Fermi 
surface; if exciton pairs within a shell of thickness A condense around the Fermi surface 
centred at kp, then the anomalous scaling (17115]) is consistent with a BCS mechanism 
with A ~ 0{fi). Everything is to be viewed in the context of an effective held theory 
valid near the QCP. In order to put this picture on hrmer footing, and also to expose the 
Fermi surface, in this paper we use Monte Carlo simulation to calculate the quasiparticle 
dispersion relation E{k), identifying kp with the location of the minumum and A with 
E{kp). The main results are summarised in Fig. Ejbelow. In Sec. Elwe present the model 
and review the main hndings of [T], then in Sec. El present the calculation of E{k). Our 
results, summarised in Sec. HI indeed support the picture of a Fermi surface disrupted by 
strong interactions, leading to the formation of a gap A increasing monotonically with 
/i. We also discuss our observation of the striking inequality < kp, and propose an 
explanation in terms of an estimate for the dynamical critical exponent z < 1. 


2 


2 Formulation and Simulation of the Model 


The model we use to describe voltage-biased bilayer graphene in terms of iVj = 4 rela¬ 
tivistic fermions is described by the following Lagrangian [1]: 


C = ('ll), 0) 


(D[V-fi] 

V 


ij \ 
D[V- -^^]) 



-h At 142. 

2g^ 


( 1 ) 


Here 0 and 0 are 4x2-component spinors each describing two Dirac flavors, with 0,0 
corresponding to electron degrees of freedom on one layer and 0, 0 holes on the other; 
V is an auxiliary field dehned on the timelike links of the lattice which approximates 
an “instantaneous” Coulomb potential governed by S-l-ld Maxwell electrodynamics; and 
j a symmetry-breaking gap parameter due to interlayer pairing. Because some weak 
interlayer hybridization is likely to be present in double-layer systems, in general j 0 
0 [5]; however we will attempt to extrapolate j —)■ 0 so that exciton condensation can 
be viewed as a spontaneous symmetry breaking U(4)®U(4)—)-U(4) |1]. The bias voltage 
is given by 2/r where fi is formally equivalent to the isospin chemical potential in QCD. 
In continuum notation the covariant derivative operator is 

D[H;/i] = h“’A + + =-D^[V;-fi], (2) 

\ u =0,...,2 J 


where a, [5 run over Nf = 2 Dirac flavors. The minimal coupling to V implies that 
00, 00 and 00 interactions are all of equal strength, which is required for the action 
to be real following integration over the fermions. This corresponds to the interlayer 
separation d —)■ 0 in the double-layer model [5]. 

In terms of staggered fermion fields living on the sites x,y of a. 2+ld cubic lattice D 
is written 


TJOxG ^(1 ^ ^ Vuxidy^x+u dy^x-v) 

v=l,2 

(3) 

where the sign factors rji^x = (—ensure a covariant weak-coupling continuum 
limit. Eq. (|3]) was also used to model electron excitations in monolayer graphene [7]. 
Note that in 2-l-ld a single staggered fermion automatically describes Nf = 2 continuum 
flavors [8]. We work in units in which both spatial and temporal at lattice spacings 
are set to unity (equivalent to setting the bare Fermi velocity vp = 1); note that for 
a non-covariant action there is no reason a priori to assume = a*, though since the 
value of the ratio at/Og is driven by UV physics it does not depend on /i and only weakly 
on j, so may be assumed constant throughout this paper. On the assumption that the 
dimension of is even, it is straightforward to show 

detA4 = det[D'^D + > 0 (4) 
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and hence there is no obstrnction to Monte Carlo simnlation nsing orthodox nnmerical 
techniqnes. 

In [T] we presented resnlts from nnmerical simnlation of the model flTlIdl) nsing a 
hybrid Monte Carlo algorithm. The conpling g~‘^a' = 0.4, where the factor a' = 
follows becanse the interaction conples charge densities, was chosen in the vicinity of the 
qnantnm critical point (QCP), althongh since Nf = A ^ Nf^ = 4.8(2) [6] it is hard to 
ascertain with conhdence on which side of the phase bonndary at /r = 0 we are sitting. 
Two principal observables were monitored as a fnnction of fi and j, namely the carrier 
density 

IdhiZ 1 - 

= - =-{i/jDaTp - (pDo(p), (5) 

and the exciton condensate 

f) In — — 

(TT) = = z(-0(/) - (6) 

Following extrapolation to the limit j —)■ 0 (so long as /rat < 0.3 at which point satn- 
ration artifacts set in and the continnnm approximation fails), both observables showed 
behavionr consistent with rising smoothly from zero as fi is increased, with 


n, oc 
(TT) oc 


(7) 

( 8 ) 


This is to be contrasted with the expectations from weak-conpling. In free held theory 
the carrier density depends on the volnme contained within the Fermi snrface; for rel¬ 
ativistic fermions kp ^ fi and hence in 2+ld n,, oc /i^. Similarly, in weak conpling we 
expect the exciton condensate to arise from electron-hole pairing with eqnal and oppo¬ 
site momenta from within a shell of thickness 2A centred on kp] hence (TT) oc A/r. 
The non-canonical scaling is taken to be a symptom of strong held hnctnations near the 
QCP. Moreover, since according to Lnttinger’s theorem ric oc even in the presence of 
interactions, we can adapt this argnment to estimate the scaling of the gap A(/r): 

(T4/(/i)) I constant weak conpling; 

A(/i) = —^-oc (9) 

^c(Ai) ^ ^ ’ near QCP. 

While the nnmerical valne for the exponent shonld probably not be taken too serionsly, 

the qnalitative diherence between the scaling predicted in weak conpling and observed 

1 

near the QCP is striking; indeed. Fig. 13 of pQ shows the ratio (4/\k)/nc extrapolated 
to j ^ 0 almost linearly proportional to fi. Since near a QCP /r is the only energy scale 
in the problem, naively A oc /r is expected. As both argnments contain assnmptions, it 
is clear that a direct calculation of A from the quasiparticle propagator (\k(0)T(a;)) is 
needed. 
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3 Quasiparticle Dispersion 


While the results outlined in the previous section are intriguing the conclusions, drawn 
from simulations of systems with both significant UV and IR artifacts on a restricted 
range of /i values, are necessarily provisional. In this paper we will present complemen¬ 
tary information through analysis of the quasiparticle dispersion relation E{k), which 
will enable identification of the Fermi momentum kp and direct estimation of the super¬ 
fluid gap A. The basic observables are the timeslice correlators in momentum space: 

C^(k,t) = ^{V.(0.0)i/:(S,t))e-‘“; = ^{<*(0,0).^(f, i))e-‘“ (10) 

X X 

where we distinguish between normal propagation of an electron or hole within a layer, 
and anomalous propagation in which an electron moving in one layer is absorbed by 
an exciton, transferring its momentum to an electron moving in the other. On a hnite 
volume in the absence of explicit symmetry breaking j = 0, the anomalous component 
Ca necessarily vanishes. 

Many momentum modes must be available in order to obtain a good resolution for 
E{k). Naively this would entail making at least one of the spatial dimensions of the 
lattice as large as possible [H]. It is much more efficient, however, to employ partially 
twisted boundary conditions [10] in which the constraint 'ipi.Lx) = with the 

angle adjustable, is implemented in the calculation of the propagator (1T0|) so that 
accessible modes have 

2Tin -|- 9x 

kx =-- 

Lx 

In practice the twist is implemented via a held redehnition so that each x-link in ((3|) is 
multiplied by a phase Treating +x and —x independently, the choice 9x = 27r/3 

permits an effective tripling of the number of available modes for a given Lx for the 
cost of an extra inversion (although equivalent statistics requires twice as many twisted 
inversions as non-twisted). For the two volumes investigated here, the accessible modes 
were kog = 0, ..., ^ . (32^) and kog = 0,^,...,^... (48^). In all cases the 

maximum accessible momentum for staggered fermions is 

Results were generated at g~^a' = 0.4 at /ia* = 0.1, 0.2 on 32^ and jaat = 0.1, 0.15, 0.2 
on 48^, with exciton source ja = 0.005, 0.01, 0.02,..., 0.05. (here a = with 

a dehned by the renormalisation prescription adopted for the non-conserved density 
■ 00 ). A hybrid Monte Carlo algorithm with 5t = 0.0025 and mean trajectory length 
r = 2 (32^) or r = 1 (48^) was used, and propagator measurements taken at the end 
of every trajectory using a randomly chosen point source. The results presented arise 
from 0(5000) measurements at ja = 0.005 up to 0(3 x 10"^) at ja = 0.05. Since the 
conventional staggered fermion mass term is absent from ([3|), the global symmetry 

^ ^ ^-iae{x)^ ^ ^-iae(x)^ ^^2) 

with e{x) = (^—iYo+xi+x 2 ^ means that C^ik, t) vanishes for t even and CA{k, t) for t odd. 
Moreover only RoCat and ImC^ survive the ensemble average. The resulting correlators 
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are then fitted to the forms m 


(13) 

(14) 


CN{k,t) = + 

CA{k,t) = 

to yield the fc-dependent amplitudes A, B, \C\ and the energy E, which is extracted 
independently from both Cn and Ca- Stable hts for CAr(f) and C'yi(f — 1) were found 
for the windows tat ^ [7,25] (32^) and tat £ [7,41] (48^). It is slightly unusual to have 
to deal with excited state contamination when htting a fermion propagator, which may 
be a symptom of the strong fluctuations near a QCP; this was found to be an issue 
particularly for k ^ kp. 



Figure 1: Amplitudes obtained from fits to (I13I14I) on 48^ at p,at = 0.1. 

It is instructive hrst to consider the htted amplitudes: Fig. [1] shows results from 
fiat = 0.1 on 48^, for two values of the exciton source j. In the normal channel the time- 
asymmetric form of Ctv is manifest, and changes in character as k increases. For small 
k the forwards propagating signal is stronger, but B/A grows with k and for kag ~ 0.2 
the backwards signal dominates. The interpretation is as follows m- for k < kp the 
dominant excitations are hole-like, and for k > kp particle-like. Increasing jd from 
0.005 to 0.05 has the effect of smearing the Fermi surface so that quasiparticles tend 
to become an admixture of both, and the disparity between A and B diminishes. The 
effect of smearing is also seen in the anomalous channel, where Ca grows steadily in 
magnitude with increasing j. In weakly-coupled models [HI E] \C{k)\ is non-monotonic 
with a maximum near kp where particle-hole mixing is strongest, but here the behaviour 
is less clear-cut. The same trends with both k and j are observed at larger /i. 
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Figure 2: Dispersion relation E{k) in the normal channel (unless stated) for various /x, for fixed ja = 

0 . 01 . 

Fig. [2] shows results from fits to the dispersion E{k) for various fi at ja = 0.01 
from fll3lll4p . The common feature is that E{k) is non-monotonic with a minimum in 
the neighbourhood where the amplitude ratio A/B ^ 1, which we have identified as 
the Fermi momentum kp. For k < kp quasiparticle excitations are hole-like, and the 
energy needed to excite them from the ground state decreases as k kp. For k > kp, 
excitations are particle-like and the opposite holds true: indeed, in this regime results 
from all three /i-values studied are plausibly consistent with being drawn from the same 
branch of the dispersion curve appropriate to the vacuum (ie. with zero bias voltage). 
Fig. [2] compares results from two volumes 32^ and 48^, and also for ^at = 0.15 with fits 
in both normal flTSll and anomlaous ffTTll channels. On the assumption that a smooth 
curve may be drawn through the admittedly noisy data, there is no evidence for any 
signihcant hnite volume artifacts, or systematic difference between the two channels. 

Fig. [3] plots E{k) for various j at fiat = 0.1 on 48^. The non-monotonicity observed 
above becomes more pronounced as j —)■ 0, and since there is little shift in the minimum 
with j, we adopt the pragmatic procedure of identifying the Fermi momentum kp with 
the value of k where E is minimum. The resulting estimates are shown in Table (H the 
quoted error is half the mode spacing on 48^, except at fiat = 0.2 where the dispersion is 
flatter and a full mode spacing is taken. The first observation is that kp is systematically 
greater than fi, consistent with the precocious saturation of ndfi) observed in [T]. This 
is discussed further in Sec. HI but already we note a discrepancy with the expectation 
kp = fi for free massless fermions with = at- Setting aside the issue of the smearing 
of the Fermi surface by an exciton gap A > 0, we can at this stage test consistency with 
Luttinger’s theorem, which states that Uc depends solely on the geometry of the Fermi 
surface characterised by kp, independent of the nature of the interactions. The third 
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Figure 3: E{k) in the normal channel for ^at = 0.1, for various j on 48^. 


fjjdt 

kpttg 


nTikp)ai 

ndfJ-J -t 0)a^ 

0.10 

0.175(22) 

0.0032 

0.011(3) 

0.0095(1) 

0.15 

0.262(22) 

0.0080 

0.023(4) 

0.0328(2) 

0.20 

0.436(44) 

0.0139 

0.068(15) 

0.0905(3) 


Table 1: Estimates for the Fermi momentum kp and comparison with Luttinger’s theorem. 


column of Table [T] shows the carrier density evaluated for free massless fermions on 
the same 48^ lattice at the reference value of /i, and the fourth with fi set equal to kp 
in the second column; the hfth column gives the value of Uc in the interacting theory 
in the limit j —)■ 0 obtained in [I]. In all cases while the values 

of 'n}^‘^^{kp) and nc(/i) are comparable, which is encouraging; however as fi increases 
'n}^‘^'^{kp) < ricifi) indicative of the difficulties in precisely locating kp due to both the 
limited momentum resolution and perhaps the absence of a sharp Fermi surface due to 
exciton condensation. 

In Fig. IHwe plot the gap A dehned as the energy E{kp). Data from all available 
hts on both volumes is shown, in both normal and anomalous channels. As before 
there is little evidence for a systematic effect with lattice volume and both channels 
yield consistent results. The dashed lines show an extrapolation to j = 0 based on the 
quadratic form 

A = Ao + aj + bf. (15) 

It is clear hmj^oA(j) ^ 0 which is direct evidence for spontaneous gap formation via 
exciton condensation for all three values of fi. The resulting extrapolations, together 
with the estimates for kp{fi) from Tabled] are shown in Fig. [S] which is the main result 
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Figure 4: Gap A vs. j for various /i. Filled symbols are from 48^, empty from 32^, and the dotted lines 

are fits to normal data on 48^ using (HSJ. 


of this paper. In the region of fi studied it is clear that A varies strongly with and is 
of the same order of magnitude as the chemical potential /r, which are notable features 
of this particular model proposed in [T], and indicative of strong interactions at the 
Fermi surface. This conclusion holds in both normal and anomalous channels, appears 
to be independent of volume, and is striking enough to be robust against uncertainties 
introduced by the ad hoc nature of our analysis (an analytic form against which to fit 
the dispersion E{k) would be very valuable), the IR and UV artifacts inherent in studies 
of lattice models with /i 7 ^ 0 [ 12 ], and lack of knowledge of the physical anisotropy at/ag. 
On the basis of the three chemical potentials studied both A and kp appear to scale 
superlinearly with /i, in qualitative agreement with (jH]). 


4 Discussion 

In this paper we have used lattice Monte Carlo simulation techniques to explore the 
quasiparticle dispersion relation in an interacting held theory with non-zero charge den¬ 
sity, and shown that for /c ~ fci? the excitations are gapped with A ~ 0(/r), and A scaling 
faster than linearly with /i. This is in sharp contrast to results from comparable studies 
of other simulable models with /r 7 ^ 0. In [9] the gap in the 3+ld Nambu Jona-Lasinio 
model (a relativistic analogue of the original BCS model) was shown to be approximately 
constant above onset, independent of and numerically much smaller than /x, consistent 
with the BCS result A ~ Af/y exp(—cA^y//i^). In QCD with gauge group SU(2) there 
is a so-called quarkyonic regime above onset where (TT) oc /x^ [ 12 ]; this is consistent 
with degenerate fermions in 3+ld with a gap A ~ 0{Aqcd) independent of /x. It is also 
very different from the result A//i ~ 0(10“^) obtained by self-consistent diagrammatic 
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Figure 5: A(j = 0) versus /r. The square symbols denote kp from Table [T] 


techniques [13], although comparable with the large values of A//i obtained in [5], where 
it was found that A depends sensitively on the treatment of screening effects, and in 
particular on the reduction of screening once the superfluid gap forms. One feature of 
our approach which does merit comparison with the treatment in [5] is that competition 
between inter- and intra-layer pairing condensates can be addressed; see Fig. 11 of |T]. 

We now return to Tabled] and the issue of why pot < kpag- In |T] it was suggested 
this is because in a strongly self-bound system the Fermi energy is necessarily less than 
the Fermi momentum in natural units. While this may be plausible for a system where 
/i A is by far the largest scale, it is difficult to see how this picture can persist in the 
regime we have been focussing on. Another possibility, which we cannot dismiss, is that 
the disparity is a lattice artifact caused by a large induced anisotropy at/as ~ 0(0.5). 
Indeed, if we assume the Fermi velocity remains close to one even in the presence of 
interactions, then the dispersion data of Fig. [2] might suggest at/ag ~ 0.3. However, 
a more compelling possibility is that E{k) is not a linear relation, but rather a power 
law characteristic of a nearby QCP. Rewrite the scaling form o as Uc oc /i^; the 
Luttinger scaling Uc oc then gives E cc k^, where guided by Fig. [2] we assume the 
relation between Ep and kp completely characterises the quasiparticle dispersion. In this 
scenario the scaling (I7|) extracted from bulk observables in [T] thus yields an estimate 
for the dynamical critical exponent 


z ~ 0.6. (16) 

Much greater numerical precision than achieved in Fig. |2] would be needed to dis¬ 
tinguish these two possibilites unambiguously. Another route would be to perform a 
“biased bilayer” study for a related 2+ld theory, the Thirring model mi. whose be- 
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haviour as a function of Nf and /i is qualitatively similar to the model here [15], but 
whose continuum action is manifestly covariant implying at = throughout. 

Finally, we note that the superlinear scaling of A(/i) is also suggestive of a power 
law A oc , with a > 1, modifying our naive expectation A oc /i. Clearly, a much more 
extensive simulation campaign is required to verify this interesting possibility. 
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